Outer boundary conditions for Einstein's field equations in harmonic coordinates 
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We analyze Einstein's vacuum field equations in generalized harmonic coordinates on a compact 
spatial domain with boundaries. We specify a class of boundary conditions which is constraint- 
I preserving and sufficiently general to include recent proposals for reducing the amount of spurious 

, reflections of gravitational radiation. In particular, our class comprises the boundary conditions re- 

' cently proposed by Kreiss and Winicour, a geometric modification thereof, the freezing- 5'o boundary 

04 ' condition and the hierarchy of absorbing boundary conditions introduced by Buchman and Sarbach. 

^ , Using the recent technique developed by Kreiss and Winicour based on an appropriate reduction 

Q ' to a pseudo-differential first order system, we prove well posedness of the resulting initial-boundary 

' value problem in the frozen coefficient approximation. In view of the theory of pseudo-differential 

, operators it is expected that the full nonlinear problem is also well posed. Furthermore, we imple- 

■ ment some of our boundary conditions numerically and study their effectiveness in a test problem 
( ' consisting of a perturbed Schwarzschild black hole. 

I— l! PACS numbers: 04.20.-q, 04.20.Ex, 04.25. Dm 

cr. 

)-<'. I. INTRODUCTION 

OXJ, 

A common way to deal with the numerical simulation of wave propagation on an infinite domain is to replace 
the latter by a finite computational domain S with artificial boundary dYj. Boundary conditions at dT, must then 
be specified such that the resulting Cauchy problem is well posed. Additionally, the artificial boundary must be as 
transparent as possible to the physical problem on the infinite domain in the sense that it does not introduce too much 
f — . . spurious reflection from the boundary surface. The construction of such absorbing boundary conditions has received 
CNJ ■ significant attention for wave problems in acoustics, electromagnetism, meteorology, and solid geophysics (see [l| for 
, a review). 

■ In this article, we construct absorbing boundary conditions for Einstein's field equations in generalized harmonic 
' coordinates. In these coordinates, one obtains a system of ten coupled quasi-linear wave equations for the ten 

components of the metric field. Therefore, ten boundary conditions must be imposed. However, it is important to 
^ I realize that not all ten components of the metric represent physical degrees of freedom, so that one cannot simply 

• ^ apply the known results from the scalar wave equation directly to the ten wave equations for the metric components. 

^\ Instead, one has to take into account the fact that the metric is subject to the harmonic condition which yields 
^ four constraints. In the absence of boundaries, it is possible to show that it is sufficient to solve these constraints 
■ - - ' along with their time derivatives on an initial Cauchy surface; the Bianchi identities and the evolution equations then 
guarantee that the constraints are satisfied everywhere and at each time. When timelike boundaries are present, four 
constraint-preserving boundary conditions need to be specified at dY, in order to insure that no constraint-violating 
modes propagate into the computational domain. This reduces the ten degrees of freedom to six. Another four degrees 
of freedom are related to the residual gauge freedom in choosing harmonic coordinates and fixing the geometry of 
the boundary surface. Therefore, one is left with two degrees of freedom which are related to the gravitational 
radiation. The challenge, then, is to specify boundary conditions at dY, which preserve the constraints, form a well 
posed initial-boundary value problem (IBVP) and minimize spurious refiections of gravitational radiation from the 
boundary. Additionally, one might want to require that gauge and constraint-violating modes propagate out of the 
domain without too much reflection. 

Constraint-preserving boundary conditions for the harmonic system have been proposed before in 0, 0, 0j [F| sjid 
tested numerically in [1, 0, 0, B H [ifl] ■ The boundary conditions of [3, Q are a combination of homogeneous Dirichlet 
and Neumann conditions, for which well posedness can be shown by standard techniques. On the other hand, the 
conditions of [5, 's'l are likely to yield large spurious reflections of gravitational radiation and probably do not give a 
good approximation to the solution on the unbounded domain. Inhomogeneous boundary conditions which allow for 
a matching to a characteristic code are also considered in [3], but in this case the well posedness of the problem is 
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not established. The boundary conditions of Q are of the Sommerfeld type, and the well posedness of the resulting 
IBVP has been shown, at least in the frozen coefficient approximation. Finally, the boundary conditions presented 
in [3] contain second derivatives of the metric fields and freeze the Weyl scalar \l/o to its initial value. As discussed 
below, this condition serves as a good first approximation to an absorbing condition. In this paper, we generalize the 
first order boundary conditions of ^5] to the full nonlinear case, and also obtain more general second order boundary 
conditions that are very similar to the conditions presented in Furthermore, we obtain a class of constraint- 
preserving higher order boundary conditions which are flexible enough to incorporate the recently proposed hierarchy 
of absorbing outer boundary conditions proposed in [l^, [l^ ■ 

There has been a considerable amount of work on constructing well posed constraint-preserving boundary conditions 
for Einstein's field equations. A well posed IBVP for Einstein's vacuum equations was presented in Ref. [3]. This 
work, which is based on a tetrad formulation, recasts the evolution equations into first order symmetric hyperbolic 
quasilinear form with maximally dissipative boundary conditions p3 , , for which (local in time) well posedness is 
guaranteed [l6j . There has been a substantial effort to obtain well posed formulations for the more commonly used 
metric formulations of gravity using similar mathematical techniques (see ^'^^ partial results). A different 

technique for showing the well posedness of the IBVP is based on the frozen coefficient principle where one freezes 
the coefficients of the evolution and boundary operators. In this way, the problem is simplified to a linear, constant 
coefficient problem on the half-space which can be solved explicitly by using a Fourier-Laplace transformation [20j . 
This method yields a simple algebraic condition (the determinant condition) which is necessary for the well posedness 
of the IBVP. Sufficient conditions for the well posedness of the frozen coefficient problem were developed by Kreiss 2lj. 
Kreiss' theorem provides a stronger form of the determinant condition whose satisfaction leads to well posedness if the 
evolution system is strictly hyperbolic. One of the key results in [2l| is the construction of a smooth symmetrizer for 
the problem for which well posedness can be shown via an energy estimate in the frequency domain. Using the theory 
of pseudo-differential operators, it is expected that the verification of Kreiss' condition also leads to well posedness 
for quasilinear problems, such as Einstein's field equations. Work based on the verification of the Kreiss condition in 
the Einstein case is given in [1, [12]. In particular, generalized harmonic gauge and second order boundary conditions 
similar to the ones considered here were analyzed in However, since in those cases the evolution system is not 
strictly hyperbolic, it is not clear if those results are sufficient for well posedness. Kreiss' theorem was generalized to 
symmetric hyperbolic systems in (23j but their treatment assumed maximally dissipative boundary conditions. On 
the other hand, the recent work by Kreiss and Winicour ^ introduces a new pseudo-differential first order reduction 
of the wave equation which leads to a strictly hyperbolic system. Using this reduction they are able to verify Kreiss' 
condition and in this way show well posedness of the IBVP for Einstein's field equations in harmonic coordinates. 

We use this frozen coefficient technique in order to analyze the well posedness of the IBVPs resulting from our 
different boundary conditions. To this end, we consider small amplitude, high frequency perturbations of a given 
smooth background solution. In this case, the problem reduces to a system of ten decoupled wave equations on a 
frozen metric background on the half space with linear boundary conditions. By performing a suitable coordinate 
transformation which leaves the half space domain invariant, one can obtain all the metric coefficients to be those of 
the flat metric with the exception of the component of the shift normal to the boundary (see also [9|). We then prove 
using the Fourier-Laplace technique and Kreiss' theorem that our frozen coefficient problem is well posed. In view 
of the existence of a smooth symmetrizer and the theory of pseudo-differential operators [24] , it is expected that one 
can show well posedness of the full nonlinear problem as well. 

This paper is organized as follows. In Sec. [Til we summarize the generalized harmonic formulation of general 
relativity. In Sec. IIIII we study the resulting wave evolution equations for the ten components of the metric on a 
manifold of the form M — [0,T] x S, where S is a three-dimensional compact manifold with smooth boundary dT,, 
and we present several possibilities for first, second and higher order boundary conditions, where here the order 
refers to the highest derivative of the metric appearing in the condition. In the first order case, in Sec. IIII Ai we 
use the harmonic constraint to impose four Dirichlet boundary conditions for the constraint propagation system. 
Next, we speciiy two boundary conditions in terms of the shear of the outgoing null congruence associated with 
the two-dimensional cross sections of the boundary surface. These conditions are a geometric modification of the 
boundary conditions presented by Kreiss and Winicour Finally, we specify four more boundary conditions with 
some absorbing properties on the gauge modes corresponding to the residual gauge freedom. In Sec. IIII 5] we present 
second order boundary conditions. One of the advantages of allowing for second derivatives of the metric is that one 
can formulate boundary conditions for the Weyl curvature scalar 4*0, which has some attractive properties. First, ^Pq 
(with respect to a suitably chosen tetrad) represents the incoming radiation at past null infinity. Second, the Weyl 
tensor from which ^'o is constructed is a gauge-invariant quantity in the weak field limit of gravity. Third, if spacetime 
is a small perturbation of a Schwarzschild black hole, as is the case near the boundary if the boundary is far enough 
from the strong field region, ^'o (with respect to a tetrad adapted to the Schwarzschild background) is invariant with 
respect to infinitesimal coordinate transformations and tetrad rotations. A boundary condition considered in the 
literature is the so-called freezing- VPo condition 0, H, [13) [H, Eli HE HI, [131 , which freezes \E'o to its initial value. An 
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estimate for the amount of spurious reflections of gravitational radiation was given in There, a new hierarchy 
Bl, L = 1,2, 3, of conditions on ^'g was also derived with the property of being perfectly absorbing for linearized 
gravitational waves with angular momentum number smaller than or equal to L. Generalizations of these conditions 
which take into account correction terms from the curvature were presented in fl3\. In order to incorporate these 
conditions into our analysis, we consider boundary conditions of arbitrarily high order in Sec. IIII Cl In Sec. IIVI we 
use the method by Kreiss and Winicour to show the well posedness of the resulting IBVPs in the frozen coefficient 
approximation. In particular, we allow for a non-trivial shift vector, which is iinportant in view of the generalization 
to the quasi-linear case. In this sense, our results generalize the work in Ref. M to non-trivial shifts and boundary 
conditions of arbitrarily high order. Next, in Sec. [Vl we obtain estimates for the amount of spurious reflections for 
the boundary conditions constructed in this article and perform numerical tests based on a perturbed Schwarzschild 
black hole. Finally, we conclude in Sec. |Vll 



II. THE FIELD EQUATIONS IN GENERALIZED HARMONIC COORDINATES 

In this section we review the formulation of Einstein's field equations in generalized harmonic coordinates [28l . [2^ 

= Dgx^ = g^'r^b , (1) 

where H'^ are given functions on M, gab is the spacetime metric and Og = — (?''''VaVf, and V^ab denote the corre- 
sponding d'Alembertian operator and Christoffel symbols, respectively. Instead of adopting the gauge condition ([1]), 
we find it convenient to choose a fixed background manifold (M, g^,j) and replace H]) by = 0, where the vector field 
is given by 

c^ = g'^'(r^6-r^a6) -ff^ (2) 

Here, H'^ is a given vector field on M and V^ab are the Christoffel symbols corresponding to the background metric g^^b- 
In the particular case where the background manifold is Minkowski spacetime and standard Cartesian coordinates are 
chosen on M, Vab vanishes, and the condition — reduces to Eq. ([1]). The advantage of using is that, unlike 
Ogx'^, it transforms as a vector field since the difference between the two Christoffel symbols, 

C^ab = r'^afc — r'^ab — '^Q'^'^ {^ahbd + Vf,/lad ~ ^dhab^ , (3) 

forms a tensor field. Here and in the following, hab = gab — gab denotes the difference between the dynamical metric 
gab and the background metric (J^f,. Since X^fJ^f, — 0, one could also replace Vchab with X^gafc in Eq- I®; however, 
we prefer to express our equations in terms of the difference field hab instead of the metric gab- A condition that is 
related to = was used in [30*1 for imposing spatial harmonic coordinates. 
The curvature tensor corresponding to the metric gab can be written as 

R'^bcd = R bcd + 2t|cC"rf]h -I- 2C"- e[cC'^ d]b : (4) 

where R bed denotes the curvature tensor with respect to the background metric g^^b- Inserting Eq. ([3]) into Q and 
using Vag'^'^ = —C'^abg'"' — C'^'abg'"^, one obtains, 

Rabcd — 2 (VjVfe/lad ~ '^d^bhac + ^d'^a^cb — '^c^ahbd^ 

+ gefC'bcC^ ad - gefC^bdC^ ac + ^ (sae^ bed Qbe^ acd (5) 

From this, we obtain the corresponding expression for the Ricci tensor, 

Rab — 



^g"^ (-VcVd/lab - VaVfc/lcd + VaVe/lM + ^b^chad) 
gefg"'^ {CacC^bd - C'^abC^cd) ~ g'"^ R cd(agb)e ■ (6) 



On the other hand, we have 



VaCh = .g^"^ (Wa^chbd - ^ 4 Vfc/lcd) - 2 C'da gbe C\f g''^ 

- gefg^'^C^bCf.d-^aHb . (7) 
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Subtracting the symmetric part of Eq. ([7]) from Eq. ([6]) we obtain 

Sab = 9'''^c^dhab-2gefg'''C\,Cfbd-^C'd(a9b)eC\fg''f 

+ 2 g'^'^R cd{a9b)e - {aHb} , (8) 

where £ab = ^'2Rab + 2V((jCf,). Using the twice contracted Bianchi identities we also obtain, 

1 



Sab - ^gabg^'^Scdj = V'VbCa + Ra'Cb . (9) 

The Cauchy problem for Einstein's vacuum equations in generalized harmonic coordinates on an infinite domain of 
the form M — [0,T] x M"^ can be formulated in the following two steps. First, specify initial data on the hypersurface 
Eq := {0} X M?. For this, let and n"" denote the future-pointing unit normals to Sq with respect to the metrics gab 
and Qabi respectively. Notice that the corresponding one-forms, and ria, are proportional to each other: Ua — a ria- 
Decompose the dynamical metric in the form 

gab = -a^naTlb + 7cd {S^a - f^^Ua) {5'^b - P'^nb) , 

where ^cd'n^ = 0, P'^ric ~ 0. The pull-back of ^ab on Eq is the metric gab induced by gab on Eqj and a and /S" are 
generalized lapse and shift. Next, solve the Hamiltonian and momentum constraints n°'{2Rab — QahQ'^'^Rr.d) = for 
the induced metric gab and the extrinsic curvature, kab- Then, solve the equation C° = which yields 0, [3l| 

a = P^Daa - a''-f''''kab ~ aiuaH" + -fabP'H''), (10) 



p = p L>bP — UbCt + a ^ 7 



Dclbd - ^Db-fcd - IbcldeH^ 



ill) 



where a dot refers to the Lie derivative with respect to n", and D denotes the connection on Sq which is induced by 
V Here, we have assumed that ViTij, = for simplicity. Equations (|10|lip can be solved by either first choosing a, 
(3°- and i?" which fixes a and /3" or by choosing a, /3, a and /3° and solving for the vector field iJ". The quantities 
a, d, /?", (3°' , "fab and kab determine the initial data for hab and hab = £^hab by taking into account that 

7q6 = {-2akcd + £pjcd) [S^a + n''na){5''-b + n'^nb). 

The second step consists in finding a solution hab on M of the nonlinear wave equation ([5]) with Eab = subject to 
the initial data specified on Sq. The results in [32, 33j show that there exists a unique solution to thisproblem, at 
least if T is small enough. Finally, we observe that the evolution equations ([5]) with Eab = imply that [J 

n-^aC" = na(2i?"'' - g^'gcdR"") + (7"'^" - n'^7''^)VcCa . (12) 

Since on Sq the Hamiltonian and momentum constraints and = are satisfied, it follows that = 0. The 
evolution system for the harmonic constraint variables, Eq. ^ with Eab — 0, then guarantees that Ca = everywhere 
on M since it has the form of a linear homogeneous wave equation for Ca with trivial initial data Ca = 0, Ca = 0. In 
the next section, we analyze the Cauchy problem when the infinite domain MP is replaced by a bounded domain E. 

III. OUTER BOUNDARY CONDITIONS 

We wish to study the evolution equations ^ with Eab = on a manifold of the form M ~ [0, T] x E, where E is a 
three-dimensional compact manifold with smooth boundary (9E. We assume that the boundary surface T — [0, T] x 9E 
is timelike and that the three-dimensional surfaces Et = {t} x E are spacelike. The cross section 5*4 — {t} x 9E 
constitutes the boundary of Ej . For the following, let n° denote the future-pointing unit normal to the time-slices Et 
and the unit outward normal to the two surface St as embedded in Ef. These vector fields are defined with respect 
to the dynamical metric gab and not the background metric gab- Therefore, 

gabu'^n'' = -1 , gabu'^s'' = , g^bS^s^ - 1 • 

The vector fields n° and s° allow us to construct a Newman-Penrose null tetrad {l"^ -, fc'^,m",m°} 

r = i= (n" + s") , fc"^ = i= {n'^ - s'^) , (13) 

^ ^iv" ^iw"), m°- ^^{v°- -iw"), (14) 
V2 V2 
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where v°- and w"" are two mutually orthogonal unit vector fields which are normal to n° and (with respect to the 
metric gab)- Notice that this null tetrad is naturally adapted to the two-surfaces St] it is unique up to rescaling of 
the real null vectors 1°" and fc" and up to a rotation i— > e^fm"-, fh°' i— > e~'^ffh°' of the complex null vectors m'^ and 
about an angle Lp. 

Since the evolution equations have the form of ten wave equations (see Eq. ([8])), we need to specify ten boundary 
conditions on T . These ten boundary conditions can be divided into constraint-preserving boundary conditions, 
boundary conditions controlling the physical radiation, and boundary conditions that control the gauge freedom. 
Constraint-preserving boundary conditions make sure that solutions with constraint-satisfying initial data satisfy the 
constraints for each < t < T. Since there are four constraints, namely = 0, and these constraints obey a set of 
wave equations on their own (see Eq. ©), there are four constraint-preserving boundary conditions. Gravitational 
radiation has two degrees of freedom, so we need to provide two boundary conditions responsible for controlling the 
physical radiation. The remaining four boundary conditions control the gauge freedom. 

In the following, we discuss several possibilities for fixing such boundary conditions. We divide them into first, 
second and higher order boundary conditions, where here the order refers to the highest number of derivatives of hab 
appearing in the boundary conditions. These families of boundary conditions are discussed next. In Sec. llVi the well 
posedness of the resulting IBVPs in the frozen coefficient approximation is proven. 



A. First order boundary conditions 

In the first order case, constraint-preserving boundary conditions are specified through 

Cc = i^ahbc - \^cha^ - i/c - 0, (15) 

where here and in the following, the notation = means an equality which holds on T only. These four boundary 
conditions are Dirichlet conditions for the constraint propagation system, Eq. ^ with £ab — 0. Since this system 
is a linear homogeneous wave equation for Ca on the curved background {M,gab), these boundary conditions imply 
(by uniqueness) that solutions of this system with trivial initial data = 0, Ca = are identically zero. Using 
g""^ = —21'^°'^'^ -\- 2m^''fh''\ the four conditions are equivalent to 

= Cc/^ = — Dimrn — DkU + D,n,nl + Dmml — Hi , (16) 
= Cck'' = — Dikk — Dkmm + D-mrnk + Dmmk — Hk , (17) 
= CcTrf = — Dikm — Dklm + Djnlk + Dmmm — H^ , (18) 

where we have defined the tensor field D^ab = ^chab E^nd where the indices /, k, m, ffi refer to contraction with fc°, 
to" and to", respectively. Notice that equation psp comprises two real- valued equations. 

Next, we consider the shear associated with the null congruence along the outgoing null vector field l'^ . This quantity 
is defined as 



'^S- (7aV-^7afc7^'')Vc/rf, (19) 



1 

2 

where jab — gab + naUb — SaSb = 2m(jjWb) is the induced metric on St- Notice that a^^^ is normal to and s" and 
trace-free; hence it has two degrees of freedom. Furthermore, it depends on first derivatives of the metric. We impose 
the boundary condition 

to^toV^',) = 92 , (20) 

where q2 is a given complex- valued function on T. In order to express this condition in terms of C^ab and background 
quantities, we first notice that the one-forms Ua and Sa are related to their corresponding background quantities ria 
and Sa by 

Ha = aria , Sa = esa + S ria , 
where a, 5 and e are functions on T with a and e being strictly positive. Next, we compute 
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Since {n^ , s^} and {ria , s^} span the same vector space, ^a'^^b'^'^c^d ='■ coincides with the second fundamental 
form of the two-surface St as embedded in the background manifold (M, g^f^) with respect to the normal vector field 



^''h. Therefore, 



.(0 



Since V2la = {a + S)ria + e Sa, is explicitly given by 



^ifc = ^7a''7f>'' (a + S)\7cnd + e VJd 



V2 



For the following, it is important to notice that while K.'fl depends on the metric fields hab, it does not depend on 
derivatives of hab- Finally, using Eq. ([3]), the boundary condition l|20p can be expressed as 



(21) 



Considering that -\/2/°9a — n'^da + s'^da, we see that the boundary conditions pB|) - P^ and (^1]) yield generalized 

Sommerfeld conditions of the form I'^VaU = q for the metric components u G {hmm , hkk , h^m , hmm} where q does not 
contain any derivatives along Z". For this reason, we specify four more Sommerfeld-like conditions on the remaining 
metric components hu, hik, him and obtain the first order boundary conditions 

(22) 
(23) 
(24) 

ni,n + 2 (g2 - ICllL) , (25) 

kll + Drnlfri + D^im — Hi , (26) 
klm ^ Djjiik ^ -^mmm i (27) 

kinifi ^ D.,jifak ^" Dfj^^jik -^k t (2S) 

where p and tt are real-valued given functions on T and qi and q2 are complex-value given functions on T, with 
q2 = m°'m!'a''ll representing the shear with respect to the outgoing null vector field P. With respect to a rotation 



Dm 




P , 


Diik 






Dllm 




qi , 


Dlnini 




2D, 


Dlrrtrh 




-D 


Dlkrn 




-D 


Dlkk 




-D 



m I— *■ e^^m of the complex null vector, we have qi 
2, respectively. 



e'^^qi and q2 i— *■ e '^q2', hence qi and q2 have spin weights 1 and 



B. Second order boundary conditions 



Next, we generalize the previous boundary conditions to second order boundary conditions, which depend on 
second derivatives of hab- The motivation for this is that such conditions can be used to reduce the amount of 
spurious refiections at the boundary surface. For example, the four boundary conditions (|15p are Dirichlet conditions 
for the constraint propagation system, Eq. ([9]) with Eab = 0, which means that constraint violations are reflected at 
the boundary [9]. Such reflections can be reduced by replacing the four Dirichlet conditions with Sommerfeld-like 
boundary conditions on the constraint variables Cc, namely 

rVaCfo = . (29) 

These conditions were used in [3] and analyzed in . As before, they imply that solutions to the constraint propagation 
system with trivial initial data vanish identically. But in contrast to the condition (|15p . the condition ([29]) allows 
most constraint violations generated inside the computational domain by numerical errors to leave the computational 
domain. "'^ In view of Eq. ([7]) these conditions yield 

— l"'!^^ aCb = —Ellmffi — Elkll + Elmfhl + Elffiml 



^ One still get reflections for plane waves with non- normal incidence, for instance. See Sec. |V] for more details. 
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— 2 C'^lClcd — C^llCcefg''^ — l°'l''VaHb , 
= l°'k^VaCb = —Eiikk — Eikmm + Eimfhk + Eifnmk 

— 2 C'^'lCkcd - C'lkCcefg"^' - I'^k'^aHb , 
= l"'w}'VaCb = —Eiikm — Eiklm + Eimlk + Eimmm 

— 2 C'^lCmcd — C^lmCcefO^^ — I'^m^'V aHb , 

where we have defined the tensor field Ecdab = "^c^dhab- 

Next, we specify the com plex We yl s calar at the boundary. Such boundary conditions have been proposed in 
the literature 0, [3, [13, [H) HI) HE EE 113 ■ particular, freezing to its initial value has been suggested as a 
good starting point for absorbing gravitational waves which propagate out of the computational domain. Recently, an 
analytic study [llj has shown that this freezing- Vl/g condition yields spurious reflections which decay as fast as {kR)~^ 
for large kR, for monochromatic radiation with wavenumber k and for an outer boundary with areal radius R. The 
Weyl scalar ^I/q is defined as 

^o = Rabcd I'^m' l^m'^. 

Using the expression ([5|) for the Riemann curvature tensor, we obtain 

2 *o = -Ell mm — Emrnll + ^Ef^i^^ 

Im ~t- 2 C ImCclni 2 C uCcmm ~\~ ^aR mini ^aR lira • 

Therefore, the six boundary conditions specified so far have the form l°-l^Wa^bU = q, where u G {hmm, hkk, hkm, hmm}, 
and q depends on zeroth, first and second derivatives of hab but only contains up to first-order derivatives with respect 
to 1°". We supplement these conditions with four similar conditions on the missing metric components hu, hik and 
him- The second order boundary conditions are then 



Eiiii 




P , 


(33) 


Eiiik 




71" , 


(34) 


Elllrn 




9i , 


(35) 


Ellmni 




Ernmll ~t" ^ E(^im)lm ~1" 2 C Im^clm 








"^C^ llCcmm ''r laR valm ~ TTT-aR Urn ~ 2'(/)o , 


(36) 


Ellkrn 




— Eiklm + Eimlk + Eiffimm ~ 2 C'^'^ I Cmcd 








ClmCcefg'^ - rm'VaHb , 


(37) 


Ellmrh 




— Elkll + Elmml + Elmml — 2 C^"^ I Clcd 








C^ll Ccefg'^ - n'VaHb , 


(38) 


Ellkk 




~ -^Ikmrfi -^Immk ~^ ^Immk ^ 2 C ; Ckcd 








C^ikCcefg"^ - rk'VaHb , 


(39) 



where p and tt are real-valued given functions on T, and qi and ipo are complex- valued given functions on T, with iJjq 
representing the Weyl scalar with respect to the Newman-Penrose null tetrad constructed in Eqs. (fT3l - fT4|) . 



(30) 
(31) 
(32) 



C. Higher order boundary conditions 

The first and second order boundary conditions constructed so far can be generalized to arbitrarily high order. Let 
L > 1 and consider the following {L + 1) order boundary conditions, 

ni''\.r^ + H'l^Va,Va,...Va,^,hcd=P, (40) 
l-H''\.r- + H'k''Va,Wa,...Wa, + ,h,d^7T, (41) 
;ai;a.^^ ;a^ + i^c^dy^^y^^ y^^^^/^^^ = q^ , (42) 

together with the four constraint-preserving boundary conditions 

rr...Z''^VaiV,,...Va^Cfc = 0, (43) 
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and the two real-valued boundary conditions 

^l^^.r^-l^m'^^TO^Va,Va,...Va^_lCcde/+6(/^cd,Va/icd,...,ValVa,...Va^/lcd;r,fc^m (44) 

where the function b depends smoothly on h^d, its derivatives of order smaller than or equal to L, and the tetrad 
vectors l"", k"" and m°. 

Boundary conditions of the form (j44p have recently been constructed in [ill , [l^ . In [ll[ , a hierarchy of boundary 
conditions Bl of this form was introduced that lead to perfect absorption for weak gravitational waves with angular 
momentum number smaller than or equal to L. This hierarchy was refined in |12l | where correction terms from 
the curvature of the spacetime were taken into account. Furthermore, as analyzed in Sec. |Vl the conditions (j43p 
should yield fewer and fewer spurious reflections of constraint violations as L is increased. The geometric meaning 
of the boundary conditions (|40l41l42p is not clear. However, their importance lies in the fact that together with the 
conditions and they yield a well posed IBVP in the limit of frozen coefficients as we shall show in the next 
section. 

We end this section by analyzing how the first order condition (|20p on the shear fits into the hierarchy (|44|) . For 
this, consider the Newman-Penrose field equation 

PVaCr — m°'WaK = — ^'o + (terms quadratic in the first derivatives of hat) ■ (45) 
Here a = m°'m!'V bla = m°'m!'fj^^l is the shear, and the spin coefficient k — m°'l^Vbla can be written as 

K = —^Dmii + (undifferentiated terms in hab) ■ 

If one could impose the boundary condition that the derivative along to" of k cancel the quadratic terms on the 
right-hand side of Eq. (|45l) . one would obtain 

*o = -rVaa, (46) 

so that the boundary condition ([20| could be thought of as the "L = member" of (|44|) . We have not found a way of 
achieving this cancellation for the general case. However, in the high-frequency limit considered in Sec. IV Al we show 
that it is possible to choose coordinates such that the condition I'^hab = is satisfied everywhere and at all times such 
that Dmii = 0. Since in the high-frequency limit the quadratic terms in Eq. (|45l) can be neglected, ([45]) then reduces 
to Eq. dMl)- 

IV. WELL POSEDNESS 

In this section, we analyze the well posedness of the IBVP resulting from the evolution equations ([5]) with £ab = 
with either the first, the second or the higher order boundary conditions discussed in the previous section. We also 
consider mixed first-order second-order boundary conditions very similar to the ones used in 0, ■ In order to do 
so, we use the frozen coefficient approximation, in which one considers small amplitude, high frequency perturbations 
of a given, smooth background solution [20l. [34|. Intuitively, this is the regime that is important for the continuous 
dependence on the data, so it is expected that if the problem is well posed in the frozen coefficient approximation, 
it is also well posed in the full nonlinear case. Using the theory of pseudo-differential operators and the symmetrizer 
construction below to estimate derivatives of arbitrary high order it should be possible to prove well posedness in the 
nonlinear case as well. 

For small amplitude, high frequency perturbations of a background solution (which we take to be gab)i the evolution 
equations ([S]) with Sab — near a given point p of the manifold M reduce to 

r\p)dcddhab = 2d(aHb) = Tab , 

where (Jabip) the constant metric tensor obtained from freezing g at the point p. Furthermore, the boundary can be 
considered to be a plane in our approximation. Therefore, the nonlinear wave equation is reduced to a linear constant 
coefficient problem on the spacetime manifold H. — (0,oo) x S, where E = {{x,y, z) e R'^ : a: > 0} is the half space. 
By performing a suitable coordinate transformation which leaves the foliation Ef = {t} x E invariant, it is possible 
to bring the constant metric Qabip) to the simple form 

g{p) = -dt^ + [dx + [3 dtf + dy^ + dz^ , (47) 

with P a constant. In order to see this, assume that g{p) is given by 

g{p) = -a^dT^ + h,j{dX' + (3'dT){dX^ + 13' dT) , 
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where a is a positive constant, a constant vector, and hijdXHX^ = N'^{dX^)'^+HAB{dX^+b'^dX^){dX^ +b^dX^) 
is a constant, positive definite three-metric. Here, A, B are equal to 2 or 3, is a positive constant, h"^ a constant 
two-vector and Hab a positive definite two-metric. Then, a suitable change of the coordinates X^ and X^ gives 
Hab — Sab- Next, the transformation = NX^, = X"^ + b^X^ leaves the domain E invariant and brings 
the three metric into the form hij = 6ij. Finally, we perform the transformation t = aT, x = Y^ , y = Y^ + 
z — Y^ + P^T, which leaves the foliation Ei invariant and brings the metric into the form ((47|) . With respect to this 
metric, the evolution equations reduce to 

[-df + 2(3dtd, + (I -(3^)01 + dl + af] = Tab ■ (48) 

For the following, we assume that the constant /3 is smaller than one in magnitude. Although this condition might 
not hold everywhere if black holes are present, it holds near the boundary since the boundary surface T is assumed 
to be time-like. 



A. First order boundary conditions 



With respect to the metric (|T7)) . we have 

n'^da = dt- (3d, 

and therefore, an adapted background null tetrad to St is 



m^da 



T2^'^ 



(1- 
(1- 

idz 



(49) 



In the frozen coefficient approximation, the first order boundary conditions 



reduce to 



(50) 
(51) 
(52) 

(53) 

(54) 
(55) 
(56) 

Notice that the derivatives along fc" in Eqs. ([51)) and (|55|) could be replaced by derivatives along the time-evolution 
vector field 5* by using k^da = [V2dt - (1 - P)l°'da] /(I + (3) and Eqs. ^ and Similarly, the term k°-dahmm in 
Eq. ([56]) can be replaced by tangential derivatives by using k'^da — [V2dt — (1 — (3)1"' da] /(I + /?) and the new version 
of Eq. ([5^ . In this way, only derivatives tangential to the boundary appear on the right-hand sides of Eqs. ([5DH5S)) . 
While this observation might be useful for numerical work it is not important for what follows. The evolution system 
(I48I50I[56|) has the form of a cascade of wave problems of the form 







P , 




l"'dahlk 












Qi , 




I'^B h 




2m''dahi.m 


+ 2 (g2 - /Cll) , 


I da hni7n 




-k'^dahu + 


m'^dahlfn + fh°-dahim - Hi , 


/a o 7 




— k°'dahijn - 


(- m^'dahik + m'^dahmrn " Hm , 


l°'dahkk 




— k°-d h - 


+ m'^dahrnk + m'^dahmk " Hk 



[-d^ + 2f3dtd, + (1 - p^)dl + dl + a; 
[a* - (1 + /3)9,] = <7« , onT, 



on Q. , 



(57) 
(58) 



where z = 1 , . . . , 10 and where the boundary data qi depends on first order derivatives of the fields u^^\ for j = 1 , . . . , i — 1 
only, at the boundary surface T. 

In the following, we obtain a priori estimates for each wave problem, Eqs. ([57)1 and (I58p . using the method in 
For this, we first remark that it is sufficient to consider the case of trivial initial data. Indeed, let u = u'^*-' be any 
smooth solution to ([57)-[58)l. Then, 



u{t,x,y,z) = u{t,x,y, z) - u{Q,x,y,z) - tdtu{Q,x,y, z) , 



t>0, (x,y,z)eE, 



(59) 
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satisfies the wave problem (|57H58|) with modified source functions T^'^^ and g^*^ and trivial initial data u(0, x, y, z) = 0, 
dtu{0,x,y, z) = for all {x,y,z) S S. Next, we show there exists a constant > such that for all 77 > and all 
smooth enough solutions w*^*' of (|57ff58|) with trivial initial data, 

vWu^Xxn + h'^'^Wluri^c, + , (60) 

where the norms above are defined as 

Jn I I , 

J e-^"* J2 \drd:^d^'d:^uit,0,y,z)\'dtdydz, 

I Q I < m 

where a = {at, axjCXy,az) G Nq is a multi-index and |a| = at + ax + ay + ctz- The important point to notice here is 
that one obtains an estimate for the norm of the first order derivatives of the solution with respect to the boundary 
surface T. Therefore, in the estimate of the i'th wave problem, the norms of the first derivatives of the fields u'^', 
j = 1, i — 1, at the boundary which appear in the norm of g*-'-' on the right-hand side of (j60p can be estimated and 
one obtains the following global estimate.^ There is a constant C > such that for all 77 > and smooth enough 
solutions hab of the initial-boundary value problem (|48l50h[551) with trivial initial data. 



J2 (^ll^-&ll'.io + l|/^a6||',i.r) < \\^-b\\lo,n + \\p\\lo.T+Mlo.T 

a,b=0 \ a,fc=0 

+ ll'Zill,%,r + lk2||',o,r + 5Ill^''lllo,r) ■ 

a=0 / 

In order to prove the estimates ([60)) . let u — it*^*' be a solution of one of the wave problems ([57H58|) with trivial 
initial data, i.e. u(0, x, y, z) — 0, dtu{0, x, y, z) = for (x, y, z) e E. Next, fix > and define 

u^[T,x,y,z) - for i < 0, (a;,7;,z) G S. ^ > 

Let X, LUy, ujz) denote the Fourier transformation of x, y, z) with respect to the directions t, y and z tangential 
to the boundary, and let {t(s, x, uoy, ujz) = Urj{£,, x, tOy, oJz), s ~ rj + i^, denote the Fourier-Laplace transformation of u. 
Then, u satisfies the ordinary differential system 

[-s^ + 2f3sdx + il~f3^)dl~Lu'^]u = F, on a; e (0, cx)), (62) 
[s- {1 + (3)dx]u^ q at a; = , (63) 



where uj = ^ujy+Lol and F and q denote the Fourier-Laplace transformations of JF^'^ and q^'^ respectively. We 
rewrite this as a first order system by introducing the variable 

i = \{dx+l^l3s)u, (64) 



where k = y^sp+tJ^ and 7 = 1/ ^Jl — /3^. With respect to this, the system can be rewritten in the form 

dxW — M{s,uo)w + f , on a; e (0,00), (65) 
L{s,uj)w ^ g , at a; = , (66) 



^ Problems which satisfy this kind of estimate together with existence of solutions are called strongly well posed in the generalized sense 
in the literature [sl. 12011 . Here, "generalized sense" refers to the fact that trivial initial data is assumed and that the norms involve a 
time integration. As illustrated above, the assumption of trivial initial data does not restrict the solution space since one can always 
satisfy it by means of a transformation of the type II59I I. provided the data is sufficiently smooth. However, since this transformation 
introduces third derivatives of the initial data into the source terms Tab, it is not clear if our results can be strengthened to obtain 
strong well posedness |20|| . which does not assume trivial initial data and where the norms do not contain a time integral. 
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where 



and 



u 

i r ' k \ F 



with s' = s/k and ui' = ui/k. Notice that |s'p + = 1. Tlie eigenvalues and corresponding eigenvectors of M are 
given by 



1 



where the square root is defined to have positive real part for Re(s') > 0. One can show ^ that in this case 
Re(\/ s'^ + 7~2w'2) > Re(s') which implies that Re{fj.-) < < Re(/x+). Therefore, the solution of l[65l - l66l) belonging 
to a trivial source term, / = 0, which decays as a; — s- cx) is given by 

■w{s,x,uj) = (Te''-''e_ , (67) 

where the constant a satisfies L{s,Lu)e-a — g, i.e. 



' + + 7-2cj'2 (j = g. (68) 



It can be shown that there is a strictly positive constant (52 > such that |s' + -^/s'-^ + 7 2^'2| > for all Re(s') > 
and all a;' e M with + = 1.* Therefore, there is a constant Ci > such that 

|w)(s,0,c^)| <Ci|g(s,c^)|, (69) 

for all Re(s) > and w £ M. According to the terminology in Q this means that the system is boundary stable. The 
key result in 0, |2l| is that this implies the existence of a symmetrizer H — H{s' ,uj'), where iJ is a complex, two by 
two Hermitian matrix such that 

(i) H{s',uj') depends smoothly on {s',uj'). 

(ii) There exists a constant ei > such that 

HM + M*H > eiRe(s)/2 , 
for all Rc(s) > and all e R, where I2 denotes the two by two identity matrix. 

(iii) There are constants £2 > and C2 > such that 

<w,Hw >> e2|wp - C'2|gP , 

for all w satisfying the boundary condition L{s,uj)w = g, where < ., . > denotes the standard scalar product on 
and I . I the corresponding norm. 



^ See Lemma 2 of Ref. Q or use the following argument: Let 7]' ,a,b be real numbers such that s' = rj' + and + 7 ^uj'^ = a+ib. 

Taking the square of the last equation yields r]'^' = ah and a* + (J'^ — rj'^ — 'y~^Lo''^)o? — rj'^S^'^ = 0, from which one concludes that 
a2>V2. 

* See the proof of Lemma 3 of Ref. or use the fact that this condition is equivalent to i''+V^£_t2_ ^ > 82 for all f G C with Rc(^) > 0. 

If If I — > 00 the left-hand side converges to 2. For finite Q this inequality follows from Lemma 3.1 of Ref. [3511 . 



12 



Using this symmetrizer, the estimate (l60l) can be obtained as follows. First, using Eq. (I65p and (ii) we have 

dx < w, Hw > = 2 < w, HdxW > 

= < w, {HM + M*H)w > +2 < w, Hf > 

> eiRe{s)\w\'-K\w\^~^\Hf\\ 
where K > 0. Integrating both sides from a; = to oo and choosing K = eiRe(s)/2, we obtain, using (iii), 

2 



Re(s) J 



wl'^dx < — 



2 

£l 



< w, Hw >\x^Q + 



eiRe(s) 



\Hf\'dx 



\Hffdx. 



Since H = H{s',uj') depends smoothly on {s',uj') and \s' 



1, there is a constant C3 > such that 



\Hf \ < C3I/I for all {s' ,u>') satisfying Re(s') > and |s'p + jcj'^ = 1. Using this and multiplying the above inequality 
by k"^ on both sides, we obtain 



00 

7? J (Ifcup + dx + {\ku\^ + \dxu\'')\^^^ < C 



j \F\'dx + \q\ 




(70) 



for some constant C > 0. The estimate (150)) follows from this after integrating over ^ = Im(s), u)y and u)z and using 
Parseval's identity. Existence of solutions follows from Eqs. ([55l[M|) and standard results on ordinary differential 
equations. 

Before we proceed to the higher order boundary conditions, we remark that the estimate (|60| can be generalized 
to the following statement. For each m = 2, 3, 4, ... there exists a constant Ci,m such that 



(71) 



for all 77 > and all smooth enough solutions m'^'^ with the property that their first m time derivatives vanish identically 
at t = 0. The latter can always be achieved by means of the transformation 



u« (t, X, y, z) = u« (i, x,y,z)-Y,- {dtfu'^^ (0, x, 0) , t > , (x, y, z) G S . 



Notice that the evolution equations then imply that the first (to — 2) time derivatives of F vanish identically at i = 0. 
In order to prove the estimate ([7T|l . we first multiply both sides of (|70p by /c^^, j = 1,2, ...to — 1. This yields the 
desired estimates for the tangential derivatives. In order to estimate the normal derivatives, we use the evolution 

equation d^u = 7^ (s^ + uj'^)u — 2f3sdxU + F and the fact that 77 = Re(s) < k and obtain 



m ra 

^Y.\k"^-^diu\''dx + Y.\^^''^^ 

n J=0 3=0 



< C„ 



m-2 

/ 77-1 \k"'''''diF\''dx 
n J=0 



m-2 



a;=0 



for some constant Cm- The estimate dTT]) then follows after integrating over ^ = Im(s), o^j^ and ui^ and using Parseval's 
identity. 



B. Second and higher order boundary conditions 



Next, we generalize the previous estimate to boundary conditions of arbitrary order m > 2. In the frozen coefficient 
limit, the evolution system with the second or higher order boundary conditions discussed in the previous section has 
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the form 

[-a^ + 2/?9A + (1-/32)92 + 92 + ^^^^^ ^72) 

[ai-(l + /3)9,]'"u« onT, (73) 

where i = 1,2, ...10 and the boundary data g'*^ depends on the m'th derivatives of the fields u'^^^ for j = 1, — 1 
only. Assuming trivial initial data, defining as in Eq. (|6ip and taking the Fourier transformation with respect 
to the tangential directions {t,y,z), one obtains the same first order system ([65]) as before, but where the boundary 
condition is replaced by 



1-/9 



with the linear operator £ = (1 — /3)s' — 7 ^fc ^i9a;. In order to rewrite this in algebraic form, we notice that by virtue 
ofEq. (pjl 

Cu\_(u\_^(Q 

where the matrix B is given by 

-2 



B = 



s —7 
-72A'' 



where A' = -^/s'^ + 7 ^^'2 }jg^g positive real part. Iterating, we obtain 



Explicitly, one finds 

"'"2V -72 A' (6^+ - 6i) + bl 

where 6± = s' ± A' arc the eigenvalues of the matrix B. Therefore, the boundary conditions can be brought into the 
form ([66]) with 

L{s, u;) = ^ (&™ + 6™, -7-2a'-1(6';^ - 6")) , 

and 

m — 1 
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1-/3 



)m ^ ttL — ± 



x=0 



The solution belonging to a trivial source term, / = 0, which decays as x — > 00 is given by 

w{s,x,u!) = cre^-"'e_ , (74) 
where the constant a satisfies L{s,uj)e-a = g. Since e_ = (1, — 72A')"'", this condition reduces to 

6> = .g. 

However, as was shown in the last section, there is a constant S2 > such that > S2 for all Re(s') > and all 
w' G R with |s'|2 + |a;'|2 — 1. Therefore, there is a constant C2 > such that 

|^i(5,0,c^)| <C2|5(s,c^)|, (75) 
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for all Re(s) > and w G M and the system is boundary stable. Therefore, the exists a smooth symmetrizer satisfying 
the conditions (i)-(iii) above and we obtain the estimate 



7/ J \w\'dx+ \W\'1^^<C 




j \f?dx+\~g\^ 





for some constant C > 0. Multiplying both sides by fc^™, using the evolution equation = 



(s + uj )u — 2(3sdxU + F and 77 — Re(s) < fc, we obtain the estimate 



oc 

/m m 



< C 



x=0 



n i=0 



m-2 
j=0 



+ kT 



(76) 



for some new constant C > 0. Using Parseval's relations and assuming that (9( w(0, x, y,z) — for all j — 0, 1, ...m we 
have 



V \Hl,n,n + \\u\\lm,T <C[V '||-F||^,™-1,0 + \\P\\lm-2,T + WqWIo.t] 

Therefore, we obtain an a priori estimate as before. 



(77) 



C. Mixed first and second order boundary conditions 



In some cases, similar estimates can be proved for combinations of first order and second order boundary conditions. 
Here, we consider the boundary conditions that are obtained by combining the first order gauge boundary conditions 
(I^^HMj) with the second order constraint-preserving boundary conditions ([3BH35]) which specify '^q. This set of 
boundary conditions was used in [1, [l^l, and we shall also use them in one of our numerical tests in Sec. IV CI As 
before, we work in the frozen coefficient approximation. 

Consider first the first order gauge boundary conditions ((50l - [52l) . Using the estimate ([71]) with m = 2, we have 



3 

E 

a=0 



{v\\hla\\l2Al + \\hla\\l2,r) < Cl 



a=0 



On the other hand, applying the estimate (|77p with to = 2 to the second order boundary conditions 
high frequency limit, we obtain 

E (vWhabWl^^n+WhabWl^^T) < C2 Y {v''\\^ab\\li,n + \\^ab\\lo,T) 

a,bG{k^m,fh} a,bE{k,m.'m} 

3 3 



a=0 



a=0 



Combining the two estimates ((75H7^ we obtain 



Y {v\\hab\\l,,n + \\hab\\l2,T) < ^^3 



a, 6=0 



Y {v-'\\^ab\\i,,n + \\^ab\\io,T) 



.,b=0 



+ Y \\Ha\\li,T + \\p\\li,T + Mli,T + Ikillli.r + ll^oll^,o,r 



a=0 



(78) 
in the 



(79) 



(80) 



for a constant C3 which is independent of 7/ > and hab- Therefore, we obtain an a priori estimate also in this case. 
Notice that we have assumed that hab and its first two time derivatives vanish at t = when deriving this result. 
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V. THE QUALITY OF THE BOUNDARY CONDITIONS 



In this section, we assess the quahty of the boundary conditions constructed in Sec. IIIII We begin in Sec. IV Al 
by computing the reflection coefficients corresponding to the different boundary conditions in the high-frequency 
approximation. Then, in Sec. IV Bl we consider a spherical outer boundary and compute the reflection coefficient for 
monochromatic hnearized waves (with not necessarily high frequency) associated with the new boundary condition 
(|20p on the shear. Finally, in Sec. IV CI we implement the first- and second order boundary conditions numerically 
and compare their performance on a simple test problem. 



A. Reflection coefficients in the liigh-frequency limit 

As shown in the previous section, in the high-frequency limit our evolution system reduces to a linear constant 
coefficient problem of the form (|72I73|) on the half-space subject to the harmonic constraint 



g"' [Vahc ~ ^Vchabj = , (81) 

where Vq = da since the background metric is constant. In order to estimate the amount of spurious gravitational 
radiation reflected off the boundary, we start with a simplifying assumption: namely, we assume that the initial data 
is chosen such that it is compatible with the harmonic constraint and that the components hu , hik and him and their 
time derivatives are zero. Notice that this is not a restriction on the physics, but rather a restriction on the choice 
of coordinates as we show next. Suppose hat is an arbitrary solution of (|72|) satisfying the harmonic constraint (|8T|) . 
Under an infinitesimal coordinate transformation x'" ~ a;" -f- parametrized by a vector field f , hab is mapped to 

= hab + 2d^a^t) ■ (82) 

In particular, h'^f^ still satisfies the harmonic constraint provided that ^a obeys the wave equation 



1. 



= ^"'VcVrfCa = 2 



. (83) 



Requiring h[i, h'^. and h[^ and their time derivatives to vanish at the initial slice yields the following conditions: 

Q = h'ii - hu + V2 [dt^i - (1 + /3)9.6] , (84) 
= h'lk = hik + [dt{S,i + ^k) - (1 + P)d,ik + (1 - P)d,^i] , (85) 

= h'im = him + ^[dtU-i'i^+P)dxU + idy+id,)^i], (86) 

= Vkh'u = Vkhu + {dl + a^) , (87) 

= Vki2h'ik -Ki) = \7k{2hik -hu) + 2V2dA^i + {dl + , (88) 

= 2Vfe/iL = 2^khim + 2V™Vfcefc + {dl + dl) U , (89) 

where we have used the tetrad fields (|49p and Eq. ([55)1 . Equations ([5 ^ - (|5S)) yield elliptic equations on each x — const. 
surface for and and can be solved provided appropriate fall-off conditions at -f ^ 00 are specified. 

Once these equations are solved, Eqs. ([MH5S)) can be solved for dt^i, dt£,k and dt^m- Therefore, it is always possible 
to choose the gauge such that the harmonic constraint is satisfied and such that initially, h'^, h'^, and h'^^ and their 
time derivatives vanish. 

The evolution equations for hab, the boundary conditions ([^ - [M)) or (pOj - H^ which, in the high-frequency limit, 
reduce to 

[dt - {I + P)dx\^"^^ u ^ , u = hii,hik,him , L>0, 

and the well posedness result derived in the previous section imply that hu = hik = him = everywhere and at all 
times. In this gauge, the harmonic constraint ([5T|) yields 



^hmrfi = 0, Mhkk = —'^khmm + 2V(mft-m)fe , Mhkm = Mnh^ 
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The first condition, together with the wave equation — Vfc^^/i^ 



implies that hmm = provided 



appropriate fall-off conditions at y 



oo are specified. Knowing h,nm. , the third equation can then be integrated 



along Z° to obtain hkm- Since 1°" is outgoing at the boundary, the initial data for hkm completely determines the 
solution. Once hkm is known, the second equation can be integrated in order to obtain h^k- 

Therefore, in the gauge where I'^hat — 0, the entire dynamics is governed by the evolution system 



[-df + 2l3dtd^ 
[dt - (1 



(1 



0, 



0, 



(90) 
(91) 



where L + 1 = 1, 2, 3, ... is the order of the boundary condition considered. In order to quantify the amount of spurious 
reflections, we consider a monochromatic plane wave with frequency w > and wave vector — (3'(— 1, tan(0), 0) 
with g > and 9 E (— 7r/2,7r/2) the angle of incidence, which is reflected off the boundary a; = 0. Therefore, the 
solution has the form 

hmm = +7 6^ i^t-p,.^) ^ (^j) ^ y, z) e S , 

where (pj) = ^(1, tan(0), 0) and 7 is an amplitude reflection coefficient. Introducing this ansatz into the wave equation 
([TO]) and the boundary condition (PTjl yields the dispersion relation 



and the reflection coefficient 



u) = q 



7 = 



/3+\/l + tan2(6') 



1 



COS! 



1 + (1 + 2/3) cos( 



L + l 



(92) 



In particular, 7 = for normal incidence and 7 —1 for modes propagating tangential to the boundary {9 ±7r/2). 
For modes with fixed incidence angle — 7r/2 < 6* < 7r/2 the square bracket is nonnegative and strictly smaller than one 
which shows that fewer and fewer reflections are present if L is increased. For /3 = and L = 0, 1, the expression for 
7 given in ([5^ agrees with the coefficients obtained in Sec. l.B of Ref. |,36i] . 

Finally, we notice that the reffection coefficient 7 depends neither on the frequency nor on the wavelength. As we 
will see in the next subsection, this is an artifact of the high-frequency approximation. Also, we would like to stress 
that the result ([92)) relies on the gauge choice l°'hab = which we adopted here; it might change for other coordinate 
choices. 



B. Reflection coefficients for the shear boundary condition 



Next, we generalize the above analysis by relaxing the high frequency assumption. Instead, we assume that close to 
the boundary surface, spacetime can be written as the Schwarzschild metric of mass M (where M denotes the ADM 
mass of the system) plus a small perturbation thereof. Assuming that the outer boundary is an approximate sphere 
with areal radius R 3> M and considering monochromatic gravitational radiation characterized by a wave number 
k ^ M~^, it has been shown in [ll[ that the freezing- boundary condition yields a reflection coefficient which 
is of the order of (fci?)"'* for quadrupolar gravitational radiation. Reflection coefficients for the higher order boundary 
conditions ([44]) were also computed in [U, [13] with the result (kR)^"^^^^^^ for waves with multipole moment £ > L. 

Here, we want to compute the reflection coefficient for our shear boundary condition (|^D)) . As shown at the end 
of Sec. nil CI it can be considered as the L — member of the hierarchy of absorbing boundary conditions under 
certain circumstances such as the frozen coefficient limit in the gauge I'^hab = 0. In view of this, one could hope for 
a reflection coefficient of the order {kR)~^ for quadrupolar radiation. Unfortunately, as we show now the reflection 
coefficient only scales as (fci?,)^^ for large kR. 

In order to analyze this, consider odd-parity linear perturbations of the Schwarzschild spacetime {M,g). Hence, 
the manifold has the form M = M x 5^ and the background metric the form 

9ab = 9ij ^2;* dx^ + (JAB dx^ dx^ , (93) 

where gij denotes a pseudo-Riemannian metric on the two-dimensional orbit manifold M, r is the areal radius and 
gAB is the standard metric on S"^ . Metric perturbations of such a background have the following form 



<5ffy = , 5gAj = Qaj , 5gAB = r'^KAB , 
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where the quantities Lij, Qaj and Kab depend on the coordinates and x^. Using this notation, one finds that to 
linear order in the perturbation, the shear (jl9p associated with a, t — const, r — const surface is given by 



Saij = , 
ScTA-j = , 



(94) 



S(JAB 



1 



P 



Kab - 2 V(^Qs)j + qab r Qdj 



where here V and V refer to the covariant derivative associated with gij and (Jab^ respectively, and Kab = Kab — 
9AB g^^Kco/'^ is the trace-free part of Kab- Using the transformation properties of the fields Lij, Qaj and Kab 
under infinitesimal coordinate transformations (see for instance Ref. [HI) one can check that SaAB is invariant with 
respect to odd-parity coordinate transformations (but not under transformations with even parity). For this reason, 
in the following, we restrict our attention to odd-parity perturbations since in this case the shear boundary condition 
(|20p has a gauge-invariant interpretation. 

Perturbations with odd parity and fixed angular momentum numbers i, m are parametrized by a scalar field k and 
a one-form h — ha dx'^ on M according to 



0. 



QAj = hj Sa 



Kab = 2kW(aS_ 



(A^B) , 



where Sa — V s F with eab the natural volume element on and Y = Y the standard spherical harmonics. 
With this notation we obtain 

SCTAB 

where /i^.*"^^ is the gauge-invariant one-form [s^ 



-PJ^-^H^aSb), 



f- 



In terms of the gauge-invariant scalar <i> obeying the Regge- Wheeler equation [3 



Mm 



6M 

T3~ 



$ = 0, 



(95) 



this gauge- invariant one- form can be computed according to /i^*"^' ^ e.y V*(r<i>), where iij is the induced volume 

element on M [stI. [sqI. |40| . Therefore, the shear boundary condition (PO)) implies the following boundary condition 
for the Regge- Wheeler equation governing the dynamics of gravitational perturbations with odd parity, 

ZJVj(r$) = 0. (96) 

For the following, we assume that the coordinates a;* = {t, r) are such that 

(^)' 

In order to quantify the amount of spurious reflectio ns g enerated by the boundary conditions (j96p . we impose these 
conditions at finite radius r = R < oo and, following |ll| . consider monochromatic quadrupolar waves of the form 



idx'dx^ = -dt^ +dr^+0 



R. 



(t,r) - ala\ (^e'" (r-t) ^ ^ ^-^k (r+t)"^ _^ 



O 



(97) 



where = ^(^r + 2/r, a\ = —dr + 1/r, fc > is a given wave number and 7 is the amplitude reflection coefficient. 
Introducing ((97)l into the boundary condition ([96]) yields (neglecting the 2M/R correction terms) 

_g2»fei?, ^3 ^ (^kR)^] {kR- i) [2i{kRf -I- 3 fci? - 3 i] = . 
Solving for 7, the amount of reflection is given by 

-1/2 



\l{kR)\ 



1 



A{kRf 



[3+(fci?)2] 



(98) 



The reflection coefficient |7(fci?)| is shown in Fig. ([T]). It can be seen from Eq. ((98)) that the coefficient decays as 
{kR)^^ for large kR. This is slower than the {kR)~^ decay we had hoped for. Therefore, it is worthwhile investing 
the effort to implement the second order boundary condition, Eq. p6p . which specifies and yields a reflection 
coefficient that decays as (fci?)"*^ when '^q is frozen to its initial value. 
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kR 

FIG. 1: Reflection coeflicient |7(A;ii)| as a function of kR for the shear boundary condition (|20p for weak, monochromatic 
quadrupolar waves with wave number k and odd parity. The reflection coefflcient is of order unity for small kR and decays as 
{kR)~^ for large kR. 

C. Numerical tests 

An ideal boundary condition would produce a solution that is identical (within the computational domain) to the 
corresponding solution on an unbounded domain. This principle was used in 10] to assess the numerical performance 
of various boundary conditions. First, a reference solution is computed on a very large computational domain. Next, 
the domain is truncated at a smaller distance where the boundary conditions are imposed. The reference domain is 
chosen large enough such that its boundary remains out of causal contact with the smaller domain for as long as we 
evolve. Finally, the solution on the smaller domain is compared with the reference solution, measuring the spurious 
reflections and constraint violations caused by the boundary conditions. 

Here we use the same test problem as in [l3| . The initial data are taken to be a Schwarzschild black hole of mass M 
in Kerr-Schild coordinates with an outgoing odd-parity quadrupolar gravitational wave perturbation (satisfying the 
full nonlinear constraint equations). The perturbation is centered about a radius rg = 5M initially and its dominant 
wavelength is A « 4M. 

These initial data are evolved on a spherical shell extending from r = 1.9M (just inside the horizon; no boundary 
conditions are needed here) out to i? = 961. 9M for the reference solution and to R = 41. 9M for the truncated domain. 
The gauge source functions Ha are chosen initially such that the time derivatives of the lapse and shift vanish, see 
Eqs. (I10lll|) . This value of Ha is then frozen in time. 

A first order formulation (in both space and time) of the generalized harmonic Einstein equations is used as described 
in Our numerical implementation employs the Caltech-Cornell Spectral Einstein Code (SpEC), which is based on 
a pseudospectral collocation method. We refer the reader to appendix A of for details on the numerical method, 
the test problem, and the various diagnostic quantities discussed below. 

Four different sets of boundary conditions are compared, 

1. the first order conditions which include the vanishing shear condition (I^S)) . 

2. the original Kreiss-Winicour 0] boundary conditions, which replace Eq. ([^5]) with 

Dl„,m = 92 , (99) 

and are otherwise identical to the previous set, 

3. the second order constraint-preserving boundary conditions with ^'o freezing, Eqs. (|33p - (|39p . 

4. the same as the previous set but with the first order gauge boundary conditions ([^^ - ([M]) instead of the second 
order ones ([55)) - ([55]) : these are the boundary conditions used in P. M. [Toj . 

Our implementation of the gauge boundary conditions differs from Eqs. (|22 p - P^ or by terms of lower 

derivative order, which were found experimentally to slightly reduce reflections from the outer boundary in the 
components l°'hab of the metric. Such non-principal terms do not affect the well posedness results of Sec. IIVI 

Fig.[2]shows the norm of the difference AU of the solution on the truncated domain with respect to the reference 
solution as a function of time. This quantity is obtained by taking a tensor norm of the differences in the metric and 
its first derivatives at each point fld\ . We normalize AU by the analogous difference of the perturbed initial data with 
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FIG. 2: Difference AW with respect to the reference solution for four different resolutions {Nr,L). Left: 1. vanishing shear 
(solid) vs. 2. Kreiss-Winicour (dotted) boundary conditions, right: 3. second order boundary conditions (solid) vs. 4. second 
order boundary conditions with first order gauge boundary conditions (dotted). 
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FIG. 3: Constraint violations C for four different resolutions {Nr,L). Left: 1. vanishing shear (solid) vs. 2. Kreiss-Winicour 
(dotted) boundary conditions, right: 3. second order boundary conditions (solid) vs. 4. second order boundary conditions with 
first order gauge boundary conditions (dotted). 



respect to the unperturbed data. The results for both versions of the first order boundary conditions are very similar. 
A first peak arises when the reflection from the outer boundary reaches the center, virherc its amplitude assumes its 
maximum because of the spherical geometry. For the second order boundary conditions, the peak is smaller by about 
two orders of magnitude. For the second order boundary conditions with first order gauge boundary conditions, 
AU appears to converge away even for the higher resolutions at late times, unlike for the first order conditions. 
Unfortunately, this is not the case for the second order gauge boundary conditions. For those, AU grows at late times 
at a rate that does not appear to depend on resolution in a monotonous way. A closer look at the data indicates that 
this growth only affects the L — 1, 2 spherical harmonic basis functions. We suspect that this is a numerical problem 
related to spectral filtering (cf. [27|); so far we have not been able to cure it. Note that AU is a gauge-dependent 
quantity because the difference norm includes the entire spacetime metric. In fact, as we shall see below, inspection 
of the errors in the constraints and in the Newman-Penrose scalar (which can be viewed as an approximation to 
the outgoing gravitational radiation) suggests that the blow-up is a pure gauge effect. 

The violations of the constraints are shown in Fig. [3) The quantity C is a tensor norm including the harmonic 
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FIG. 4: Difference of ^'4 with respect to the reference solution for two different resolutions. Left: 1. vanishing shear (solid) 
vs. 2. Kreiss-Winicour (dotted) boundary conditions, right: 3. second order boundary conditions (solid) vs. 4. second order 
boundary conditions with first order gauge boundary conditions (dotted). 

constraints as well as the additional constraints arising from the first order reduction of We normalize C by 
the second derivatives of the metric so that C ^ 1 means that the constraints are not satisfied at all. The constraint 
violations converge away with increasing resolution for all the boundary conditions. This is what we expect because 
all the boundary conditions we considered are constraint-preserving. 

One of the main objectives of numerical relativity is the computation of the gravitational radiation emitted by a 
compact source. Hence it is important to evaluate how the boundary conditions affect the accuracy of the extracted 
waveform. To this end, we compute the Newman-Penrose scalar ^4 on an extraction sphere close to the outer 
boundary (at i?ox = 40M)'5. The tetrad we use agrees with the one given in Eqs. (|13p - (|14p when evaluated for the 
background spacetime (see [To'] for details). Strictly speaking, ^'4 only has a gauge- invariant meaning in the limit 
as future null infinity is approached but since our computational domain does not extend to infinity we can only 
evaluate 4*4 at a finite radius. However, 5*4 is gauge-invariant with respect to infinitesimal coordinate transformations 
and tetrad rotations on a Schwarzschild background, so errors in ^'4 due to gauge ambiguities should be very small. 
Fig. m shows the difference of ^'4 with respect to the same quantity obtained from the reference solution at the same 
location. We normalize |A^'4| by the maximum in time of |\E'4| at the extraction radius. Again, both versions of 
the first order boundary conditions show very similar numerical performance. Clearly visible is a first peak arising 
when the outgoing wave passes through the extraction sphere. Some of it is reflected off the boundary and excites the 
black hole, which then emits quasinormal mode radiation of exponentially decaying amplitude-a feature also visible 
in Fig. m The reflections are much smaller for (both versions of) the second order boundary conditions (about an 
order of magnitude at the first peak and two-three orders of magnitude later on). Unlike for the first order conditions, 
their |A4'4| decreases with increasing resolution, at least at late times. 

Finally we estimate the reflection coefficients for the various boundary conditions numerically and compare with 
the analytical predictions. As a consequence of the results of Ref. [ll|, the reflection coefRcient can be approximated 
by forming the ratio of the Ncwman-Penrose scalars ^0 and ^"4 at the outer boundary, 

\l{kR)\^^^^+0{kR)-\ (100) 

where k is the wavenumber and R is the boundary radius. For the vanishing shear boundary conditions (|25p . we 
found in Sec. IV Bl 

|7(A:i?)| = i(fci?)-i +0(fci?)-2, (101) 



^ We decompose the Newman-Penrose scalars with respect to spin-weighted spherical harmonics on the extraction sphere and only display 
the (by far) dominant mode [id . 
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FIG. 5: Comparison of the time Fourier transform of the measured 'l'o(i) (extracted 1.9M in from the outer boundary) 
with the predicted value using the reflection coefficients derived in Sec. IV Bl and in [ll]]. Left: 1. vanishing shear (solid) 
vs. 2. Kreiss-Winicour (dotted) boundary conditions, right: 3. second order boundary conditions (solid) vs. 4. second order 
boundary conditions with first order gauge boundary conditions (dotted). Top: outer boundary radius R — 41. 9A'/, bottom: 
R = 121. 9M. 

whereas for the freezing- condition (|36p . we have the much smaller reflection coefficient [ll| 

\j{kR)\^ ^{kR)-* + 0{kRy^. (102) 

In Fig.[5l we compare the measured ^I'o with the predicted value obtained from Eq. (jlOOp . using the measured and 
the above analytical expressions for the reflection coefficients. A Fourier transform in time has been taken in order to 
obtain plots vs. wavenumber k. The agreement is rather good, roughly at the expected level of accuracy 0{kR)~^. 
The leveling off of the numerical ^'o for large k is likely to be caused by numerical roundoff error (note the magnitude 
of \E'o at large k). The plots also indicate that the reflection coefficients are virtually the same for both versions of 
the second-order boundary conditions (as expected since they only differ in the gauge boundary conditions) , and that 
the reflection coefficient of the original Kreiss-Winicour boundary conditions agrees with that of our vanishing shear 
conditions. 

Summarizing, both versions of the first order conditions (those including the vanishing shear condition (j25p and 
the original Kreiss-Winicour conditions) performed very similarly in our numerical test. In contrast, the second order 
conditions caused substantially less spurious reflections from the outer boundary. 
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VI. CONCLUSIONS 

In this paper, we have derived various sets of absorbing and constraint-preserving boundary conditions for the 
Einstein equations in the generahzed harmonic gauge. We divided them into first, second and higher order boundary 
conditions where the order refers to the highest number of derivatives of the metric fields appearing in the boundary 
conditions. The first order boundary conditions are a generahzation of the conditions considered by Kreiss and 
Winicour [Bj and specify the shear of the outgoing null congruence associated with the two-dimensional cross sections 
of the boundary surface. Our second order conditions enable one to fix the Weyl scalar 'i'o at the boundary. Although 
there is a gauge ambiguity in the definition of "^/q at finite radius, these conditions allow, in some sense, control of the 
incoming gravitational radiation. This is important for simulations aimed at the far-field extraction of gravitational 
waves emitted from compact astrophysical sources. Furthermore, we could for example study the critical collapse of 
gravitational waves by starting with Minkowski spacetime and injecting pulses of gravitational radiation through the 
outer boundary with different amplitudes [3, [2^ . Finally, we have considered higher order boundary conditions which 
comprise the hierarchy of absorbing boundary conditions Bl and Cl discussed in pTi . [l^ . As was shown in these 
references, Bl and Cl yield fewer and fewer spurious reflections of gravitational radiation as L is increased. 

In Sec. IIVI we have analyzed the well posedness of the IBVPs resulting from our different boundary conditions. 
In order to do so, we considered high-frequency perturbations of a given smooth background solution in which case 
the problem reduces to a system of ten decoupled wave equations with boundary conditions on a frozen background 
spacetime. By means of a suitable coordinate transformation, we have reduced the background metric to the flat 
metric, with the exception of the component of the shift normal to the boundary. Using the technique of Kreiss and 
Winicour [s^] which is based on a reduction to a pseudo-differential first order system and the construction of smooth 
symmetrizer, we then have shown that the resulting IBVPs are well posed in the high-frequency limit. In view of the 
theory of pseudo-differential operators [13] and the fact that we obtain estimates for derivatives of arbirtrary order 
it is expected that the full nonlinear problem is well posed as well. Our results thus generalize the work of Ref. M 
to non-trivial shifts and boundary conditions of arbitrarily high order. They also strengthen the result of Ref. [9|, 
where boundary stability but not well posedness was proved for a first order version of the generalized harmonic 
Einstein equations derived in We remark that our results imply well posedness of such first order formulations 
provided that the evolution system of the additional constraints related to the first order reduction (supplemented 
with suitable constraint-preserving boundary conditions) is well posed. For a recent proof of well posedness for the first 
order boundary conditions which is based on integration by parts, and which does not require the pseudo-differential 
calculus, see [41|. 

In order to study the quality of the different boundary conditions considered in this paper, we have computed 
the amount of spurious gravitational radiation reflected off the boundary in the high-frequency approximation in 
Sec. |Vl We have shown that fewer and fewer reflections are present if the order of the boundary conditions is 
increased. In addition, we have generalized that analysis without the high frequency approximation for odd-parity 
linear gravitational waves with wavenumber k propagating on the asymptotic region of a Schwarzschild background. 
For the case of a spherical outer boundary of areal radius R with the shear boundary condition the reflection coefficient 
has been found to scale only as (kR)~^ for large kR which is much slower than the {kR)~'^ decay calculated for the 
freezing- boundary condition . Finally, we have performed numerical tests of some of our boundary conditions 
similar to the ones presented in [lO|. The initial data were taken to be a Schwarzschild black hole with an outgoing 
odd-parity quadrupolar gravitational wave perturbation. The first order boundary conditions (with our modified 
vanishing-shear condition) performed very similarly to the original conditions considered in Q. In contrast, as 
expected from the analytic considerations, the second order conditions caused substantially less spurious reflections 
from the outer boundary. A numerical implementation of the higher order boundary conditions is beyond the scope 
of this article and will be presented in future work. 
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